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ABSTRACT 

The slowing down of fast neutrons was analyzed by a 
multi-group method of discrete time and energy states 
coupled wdth a spatial harmonic expansion method to deter- 
mine the neutron density in a homogeneous, isotropically 
scattering slab. Five neutron source geometries were 
studied for both a fissioning and a non- fissioning system. 

Numerical results were obtained for the neutron flux, 
mean neutron energy and the neutron spectra for the one 
dimensional system using a harmonic mode expansion of up 
to six terms to determine the time -energy-space dependence. 
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I. INTRODUCTION 



Existing and future experimental work in the study of 
pulsed fast-neutron assemblies require the ability to de- 
termine the time dependent solution of the neutron trans- 
port equation accurately in the time period extending from 
the nanosecond to the millisecond range. This work allows 
the experimenter to compare both measured response and the 
analytical predictions for the overall improvement of each. 
Safe operation and control of fast nuclear reactors demands 
a fast, accurate mathematical simulation. 

Current research on the response of fast assemblies 
to pulsed neutron sources considers that the time depen- 
dence can be represented by an exponential decay of the 
form: 

T (t) = ^ exp (-A^t) 
i 

where the constants and can be related to the compo- 
sition and geometry of the system. A common assumption is 
that all terms above the first or fundamental will decay 
very quickly and that all measurable quantities can be 
described by a single fundamental or pseudo - fundamental 
value for A^. Experimental effort is being directed to an- 
swering the question of whether or not such a pseudomode 
does in fact exist, and if not, to understand why not. 

This work will investigate several of the aspects of 
this problem. 
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In their movement through a material, neutrons under- 
go interactions, scattering or absorption, with the atoms 
and suffer changes in direction and velocity as the result 
of both elastic and inelastic collisions. IVhile this is a 
random process for an individual neutron, the overall pro- 
cess for considering the time history and motion of the 
bulk of the neutrons is ameniable to a statistical analy- 
sis. The Boltzman time -dependent transport equation: 

1^ (r,v,n,t) = -vfi*VN(r,v,fi,t) - vE^ (v) N (r , v ,fi , t) 



+ 



J ^v' E (v* ) N(r ,v ' ' , t) • 

^ V 



f (v' ,n'->v,fi)da’dv’ + SCr,v,n,t) 

can describe the net processes occurring in any elemental 
volume of a physical system to account for all changes in 
the density of neutrons, subject to the usual conventions: 



r 



Position of the volume element dV. 



V 



Speed . 



Direction. (The vector velocity would 
be defined as v = vft) . 

N(r,v',fl',t) Density of neutrons at time t having 

velocity v ' ft ' . 



E(v') 



Speed dependent cross section for neu- 
trons with speed (v ' ) 



f (v ' ,ft ' ->v ,ft) Probability that a neutron having vel- 

ocity (v'ft') will interact to become 
one velocity (vft) . 
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The neutron density variation v:ith time considers: 

(a) leakage o£ neutrons: 



-vfi • VN(£,v,fi,t) , 



(b) losses due to interactions (scattering or absorp- 
tion) in the volume element: 



(c) scattering of neutrons into the element phase 
space (posit ion-velocity-direct ion) due to interaction of 
neutrons having other initial parameters: 



neutrons introduced into the element of phase space from a 
source in that element: 



Complete derivations of the neutron transport equation 
can be found in Tait [1] and Davison [2] and is discussed 
in detail in many other texts on nuclear reactor physics. 

During the past twenty years, numerous methods and ap- 
proximations have been formulated to solve this intergro- 
differential equation, either in a closed analytical form, 
or by a numerical approximation, for both the time dependent 
and time - independent cases. l^ile a major portion of this 
effort has been in obtaining the solutions applicable to 
thermal reactors, uranium-235 systems, the need to conserve 



-vE^(v) N(r ,v,fi,t) , 




S (r ,v , t) . 
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fuels into the distant future requires fast breeder reac- 
tors. Adequate understanding of the complete time, space 
and energy dependence of fast neutron populations requires 
the development of new analytic and simulation models for 
the fast neutron transport and fission problem. 

Analytical solutions are often based on extension of 
the existing models for thermal systems. These solutions 
usually assume 1/v or constant neutron cross sections, 
with the age-diffusion theory [3] considering the steady 
state, time independent cases and with perturbation analy- 
sis for time dependent solutions [4] . These provide repre- 
sentations of the steady state systems, and simple kinetic 
models for fast reactors. 

(1) Numerical methods generally deal with the time 
independent solutions of the transport equation by two prin- 
cipal methods: 

Multi-group theory accounts for all energy dependent 
parameters by dividing the range of neutron energies into 
some finite number of intervals and defining average values 
of the energy dependent terms for each energy group. This 
method is frequently used to study the energy and time re- 
lations and the energy and position dependence [5]. This 
method lias been used to consider up to several hundred groups 
between fission energies and the thermal condition [6]. 

(2) . Monte Carlo is a statistical method to consider 
the interactions of individual neutrons to determine prin- 
cipal integral or system values for several spatial regions 
using multiple energy groups. The study of the individual 
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neutron "life history” from its birth to its eventual loss, 
allows determination of such parameters as the effective 
multiplication factor (Kg££) . Computational requirements 
for this method limit its usefulness mainly to time inde- 
pendent calculations, and to a limited extent [7], the time 
dependent problem using very few velocity groups. 

The description of the slowing down and thermalization 
of neutrons as a probabilistic Markov chain process was 
proposed by Perkel [8] in 1960, and expanded to consider 
the thermal energy region (0.00002 to 1.034 eV) by Ohanian 
and Daitch [9] for the study of the time dependent neutron 
spectra. Jenkins and Daitch [10] have extended this tech- 
nique to formulations of the study of a pseudo-fundamental 
mode decay for time-energy response of pulsed fast systems, 
assuming that the system response can be represented by 
the lowest Fourier spatial mode of the time dependent dif- 
fusion equation for two energy regions: .001 to 1.0 eV 

and 1.0 to 1000 Kev. 

Williamson and Albrecht [11] have extended the Markov 
chain method of treatment for the slowing down problem to 
cover the continuous range from 10.5 MeV to thermal energies. 

Menzel, et at, [12] have employed the method of Ohanian 
and Daitch to provide a complete harmonic expansion method 
for the study of the space- time-energy response of pulsed 
systems .in the interval 0.01 to 1.0 eV. 

The purpose of this thesis is to apply the numerical 
methods of the Markov chain process to describe the energy 
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and time response of pulsed fast systems, with the harmonic 
expansion techniques for the spatial dependence similar to 
Menzel, in the interval 10.5 Mev to thermal energies. 

Current work on the pulsed neutron response of fast 
assemblies, considers that the time decay can be repre- 
sented by some exponential function of the form: 



where the constants and can be related to the compo- 
sition and geometry of the assembly. 

In this study, the fundamental and higher harmonic 
space modes are combined with energy-time expressions to 
determine the effectiveness of this technique is simulating 
the experiments conducted in the laboratory. Experiments 
suggest that the decay of the neutron population is often 
not best described by the pseudo -fundamental mode. Hope- 
fully, this work will aid in determining if this is due 
to the effects of a constantly varying velocity spectrum, 
or to the persistance of higher spatial modes of the time 
response remaining for longer periods of time. 



T(t) 




1 



13 




# 



II . THE SLOWING DOWN EQUATION 



A. THE BOLTZMAN TRANSPORT EQUATION 

Starting with the general form o£ the time dependent 
transport equation, reasonable assumptions \>:hich will per- 
mit this intergro -differential equation to be restated as 
the slowing down equation are: 

(1) The neutron density N(r,v,fi,t) is isotropically 
distributed in at all points in the physical media. 

(2) All sources and scattering kernels are isotropic. 

(3) No delayed neutron sources. 

(4) The physical media is composed of locally homo- 
geneous isotropic materials. 

(5) Spatial dependence is included within the bounds 
of the diffusion approximation. 

(6) Boundry conditions are that the neutron density 
is zero at the extrapolated boundry and that the extrapo- 
lation distance is constant for all energies. 

Conditions 1, 2 and 4 permit the equivalence of the 
expression for the local streaming or leakage of neutrons 
from a volume element to be described as: 

- vfi • VN(r , V t) = D(v) V^N(r , V t) = -B^D(v) N(r^,v,fi,t) 

where B^ is called the geometric buckling. 

Within this analysis, a source that is a delta function 
in time will be considered as the mathematical equation of 
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the pulsed neutron source with a very small time duration, 
which can be written as: 



S(r,v,^,t) = S ' (r ,v , 6 (t) ) /4 tt . 

Adjusting all constant terms to account for the isotropic 
conditions of the geometry under study, and to match the 
primary condition that the neutron density at a time (t=0) 
will be equivalent to the value of the source strength at 
that time we have: 

N(r,v,0) = S(r,v,0) = ^ 

n 

where solutions of the form N(r,v,t) serve as Green's func- 
tions for the pulsed conditions with an arbitrary time and 
spatial parameter. 

Starting with the basic assumption that the neutron 
density at any time can be described as an infinite sum of 
solutions of the form: 



N(r,v,t) = F^(v,t), 

n 



it is possible to transform each of the terms of the general 
transport equation to the slowing down equation in the fol- 
lowing manner: 



nr M(r,v,Q,t) 
^ 3t 






3F^(v,t) 

3F“ 
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(2) -vfi* VN(r = - D(v) V^N(r , V , t) 



= -D(v) Y, 
n 



(3) Jv'E^(v') jN(r,v' ,t)C^,£(v’ = 

J~v' L^(v'^v) Yj F(v',t)dv* 



n 



which reduces to the now modified form of the slowing down 
equation : 



y R.Cr) 

4—1 — 

n 



— F 

4- 

O C 



n 



(v,t) + S^(v) 







V ' ( V ' ->-v) • 




F^(v,t)V^R^(r) = 



S'(r,v,t) 



substituting for the spatial -velocity dependent delta func- 
tion source, the eigenfunctional expression of the slowing 
down equation is obtained: 




+ ^t(v) Fj^(v,t) 




V * E^ (v ' -^v) • 



Fj^(v' ,t)dv' 



D(v) 



F^(v,t)V^R^Cr) 



ZRnd) Fn(v.O). 
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Consolidating all expansion functions of the index-n in 
the single summation, 



I 



n 



R„Ci) 



9F Cv,t) 

at 



+ +J*v'E^(v’ v) 



Fj^(v’,t)dv’ - Fj^(v,0) ^ + D(v) Fj^(v,t)V2R^(r) 



= 0 



it is now possible to separate the spatial dependence from 
the energy time relations via the harmonic buckling factor, 
, to two expressions: 



V^R (r) 

■ Rn(l)” " ^ 



and : 



-D(v)Bj F^(v,t) = ^ Fj^(v,t) + 2^(v) ^j^(v,t) 

+ J v'Eg(v'-v) F^(v,t)dvdt - Fj^(v,0) 

Final clearing and rearrangement of all terms yield 
the decoupled expressions: 



V^R (r) - B^ R (r) = 0 






v' 



.+ / VS (v'->v) F (v',t)dv' = 0 

S IT 



with Rj^(x) ^ time independent function and F^(v,t) a space 
independent function. 
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B. MOD-5: THE DISCRETE STATE APPROACH 

In the separated form of the slowing down equation, 
the velocity- time harmonic term must satisfy the same form 
of the intergro-dif ferential equation as the space inde- 
pendent neutron density function. 



Any method that provides a solution to the space independent 
slowing down equation can be applied to solve the velocity- 
time harmonic equation. 

Williamson and Albrecht [11,13] have developed a sto- 
chastic model for neutron moderation that provides a numer- 
ical solution of the slowing down of fast neutrons in a 
finite media. Since the slowing down process is a continu- 
ous time process in which collisions may occur at any time, 
having a continuous range of energy transition possibilities, 
it may be classified as a continuous time, continuous state 
Markov process. If the neutron cross sections remain con- 
stant with time, the slowing down may be classed as a sta- 
tionary Markov system. 

Application of the formalism of the Markov process cal- 
culations provides the capability to perform straightforward 
computer- solutions for a discrete time-velocity model of the 
space independent equation. The probability that a neutron 
will be in a finite width velocity state, v^ , at a time, t^ , 



^ N(v,t) = -vD(v) B^N(v,t) -vE^ (v)N(v,t) 




N(v' ,T) dv' + S(v,t) . 
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is calculated by a computer program, MOD-5, which follows 
the slowing down of neutrons by determining a probability 



density parameter F(v^ 

neutron will be in the 

time t.=nAt. 

J 



tj) which is the probability that 
velocity state at the discrete 



a 



The Markov process describes the variation of a prob- 
ability density vector, s(nAt), which describes the proba- 
bility that a neutron will be located at some state in the 
system during its entire lifetime from birth to death. 

This state vector is given as ; 



s(]nAt) ^ 3 » Sij, ••• , 

where s^ is the probability that the neutron is in the state 
bounded be velocities v^ and v^ (where v.<v^_-,) at the 
time nAt. This method presumes that the initial state vec- 
tor, neutron source, FCO) , can be described for time t=0. 

The problem description necessary for application of a 
MOD-5 type of approximation are: 

(a) The velocity range of interest is divided into "M” 
discrete states. 

(b) Collision physics is applied to construct an array, 
P(i,j), of the discrete Markov transition probabilities of a 
neutron going from velocity state-i to state-j in the dis- 
crete time period At. 

I'fhile the details of the Markov transition matrix are 
discussed in several sources [14] , the main considerations 
are that the array for the slowing down approximation is 
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upper triangular with only positive elements, with the 
elements of a given row summing to unity, which states that 
all possible outcomes of the process are considered. 

If the initial Markov transition matrix describes the 
probability for a transition within a specified time inter- 
val, At, the probabilities for a transition at a time {nAt} 
can be determined as : 

Vt = (PAt’"- 

The probability of a neutron being in any state at a time 
(-nAt) can be determined from the (n-l)th step in the ma- 
trix product as: 

s(nAt) = s((n-l)At) ♦ PAt. 

Both the program MOD- 5 and the harmonic expansion 
method employ several series of notations and can result in 
some conflict with existing literature. In this work, the 
following definitions and terms will be used: 



s (nAt) 


the probability density vector of a 
neutron being in any state in the sys- 
tem at the time nAt. 


s. . 
1 ,1 


the probability that a neutron will be 


in a state "i" at a time t ^ , 


' the time 


and total density vector can be written as: 


t. - 
3 


= nAt 


s{tp 


’^2,j ’ • 
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While this form of notation applies to the case of a 
single decay mode, the expansion to a multiple harmonic 
expansion method requires the addition of another indexing 
parameter-n to indicate the particular harmonic mode. This 
results in the probability density expressions being written 
as : 



n,i ,3 






probability of the neutron fraction 
of the n-th harmonic mode being in a 
velocity state v. at the time t.. 

n-th harmonic term of the time veloc- 
ity dependence of the neutron proba- 
bility density. 



The contribution of each harmonic term of the neutron 
density is the product of the spatial component and the 
probability density component: 



N (x,v.,t.) = R (x) F (v.,t.) 



where the space component, R^(x) , is a measure of the frac- 
tional contribution to the neutron population of the n-th 
harmonic mode, and the term, is the probability 

of a neutron being in the velocity state v^ at time t^ for 
the n-th mode of the time decay harmonic. 

In the following discussions, the equivalence of the 
terms is described as: 



= F 



(v. , 
n ^ 1 ’ 



t . ) Av. 

r 1 



where Av^ is the width of the velocity state about the value 

V- . 

1 
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III. GEOMETRICAL MODELS 



A. SOURCE FUNCTIONS 

The main objective o£ this thesis is the investigation 
of the harmonic mode expansion method for the space, time 
and energy dependence of the neutron density in a physical 
system. Simple source functions representative of actual 
systems are desired to demonstrate this method. Three 
categories of source type were considered: a finite volume 

source and a first collision source as interior sources in 
the system, and a wide beam exterior source condition. 

The test geometries are: 

(a) Finite Volume Source - a finite volume element, 
located on an axis of symmetry and at an off axis position 
is considered to represent a hole or cavity in the slab. 

This might be envisioned as a beam port or localized neu- 
tron generator capable of providing a pulse of monoener- 
getic neutrons in a small volume in the slab. This is 
represented as a fourier square pulse on the x-axis of the 
system. 

(b) First Collision Source - this exponential dis- 
tribution of neutrons in the system is considered as a more 
realistic approach for the response of a physical system 
than the fixed source in volume condition. Rather than as- 
sume that the initial burst of neutrons would remain in the 
fixed volume element, it can be argued that they would ini- 
tially travel throughout the system with an exponential form. 
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As with the finite volume type source, this source type is 
studied for both the symmetric and unsymmetric positions. 

(c) Exterior Source Beam - a higher level of com- 
plexity in the source geometry selection process is the 
physical situation where the assembly under study is sub- 
ject to a wide beam of fast neutrons from an external 
point. This can be considered as a very simple approxima- 
tion to a first collision response from a plane source of 
neutrons on the surface of the slab. IVhile this may be 
considered to follow an exponential type distribution 
through the slab, it is approximated as a straight line or 
ramp function within the confines of the slab. 

B. CHOICE OF DIMENSIONALITY 

The principal goal of this work is the development of 
a useful mathematical model to determine the spatial depen- 
dent response of a physical system to a burst of fast neu- 
trons. In this determination, it will be necessary to refer 
to the neutron density in terms of three main parameters - 
time, space and energy. The derivation of all equations 
and principal relationships will be written for the time, 
position and velocity of the neutrons. Equivalent expres- 
sions can easily be obtained using either the neutron ener- 
gy (E) or lethargy (u) , using the relationships: 




u = -ln(Eo/E) = -2 ln(vo/v) . 



23 



CENTRAL AXIS SQUARE PULSE 



. FIGURE 1 



OFF AXIS SQUARE PULSE 
SOURCE GEOMETRY SQUARE PULSE 
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CENTRAL AXIS FIRST COLLISION -EXPONENTIAL SOURCE 
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FIGURE 2 SOURCE GEOMETRY FIRST COLLISION SOURCE 
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FIGURE 3? SOURCE GEOMETRY EXTERIOR SOURCE RAMP FUNCTION 
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Employing the appropriate expressions, it can be 
easily shov/n that equivalent terms for the neutron den- 
sity, N(r,v,t), the time -space -velocity form can also be 
modified to give the time- space -energy or time-space-le- 
thargy form. While the velocity dependent form is used 
for the purposes of demonstration within this paper, all 
calculations performed by the programs MOD-5 and MIL-6 em- 
ploy the lethargy dependent forms. 

In the determination of numerical solutions for the 
slowing down equation, reasonable methods and simple sys- 
tems are desired to demonstrate the harmonic expansion 
method for the time, velocity and spatial expressions, 
where the method is general 'enough to be easily revised to 
handle more difficult and complex problems. As shown pre- 
viously, the neutron density can be determined by the sum- 
mation of a harmonic series whose terms are the product of 
a space function and a time -velocity function. 

N(r,v,t) = ^ 

n 

where the spatial dependent function, R(£) , satisfies the 
differential equation: 

V^R (r) - B^RCr) = 0. 

In order to retain the flexibility to study and compare the 
response of physical systems to various source types, the 
selection of the dimensionality for problem selection posed 
some difficult choices. 
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The spatial functions must form an orthogonal set, a 
condition easily met by the Fourier expansion methods for 
the one dimensional problem using trigonometric functions 
for the cartesian geometry, but requires the use of regular 
Bessel functions for cylindrical geometry, and the spheri- 
cal Bessel and spherical harmonic functions for the spheri- 
cal geometry. A practical consideration was to obtain a 
balance demonstrating the harmonic expansion application 
with reasonable computing requirements. The more complex 
the spatial functions, the longer computing time and large 
core storage space required. The trigonometric functions 
can be determined quickly, while the existing methods avail 
able for the Bessel function calculations require approxi- 
mately three times as much computing time. 

In all three of the principal geometrical coordinate 
systems (cartesian, cylindrical and spherical) , studies of 
sources located at an axis of symmetry resolve themselves 
into some form of the one dimensional problem. Hoi^rever, 
the study of sources located off the axis of symmetry of 
the system, while seeming reasonable conditions to analyze, 
can lead to some very messy mathematics, such that the util 
ity of the method can become lost in the details. While 
these problems can be considered, the one dimensional case 
would allow easy comparison of all symmetrical and non- 
symmetrical cases without initially confusing the reader or 
user of. this technique prior to demonstrating the method. 
After accomplishing that aim, the multi-dimensional cases 
(two or three dimensions) can be investigated. 
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A practical matter in the consideration of the har- 
monic expansion method involves the retention of data 
within the computing machinery while running the basic 
problem or one of the many variations that can be consid- 
ered by this method. 

An important consideration while developing the com- 
puter programming for these calculations involved the se- 
lection of methods to transfer data between the two 
computational routines. The program MIL-6 was written to 
process data input from punch cards, nine-track magnetic 
tape, or the IBM-2314 series of magnetic disk. The pro- 
gram MOD-5 was run sequentially to determine the velocity 
and time characteristics of the physical system for each 
harmonic mode of the buckling constant (B^) • This gener- 
ally requires one minute of computing time for the IBM- 
360/67 with 70 velocity states, requiring 200k bytes of 
core storage for each harmonic mode. Thus, a six-mode ex- 
pansion would require approximately six minutes of central 
processor time. 

Assuming that the problem under study is going to 
deal with 150 discrete velocity states , 20 time step values 
and a minimum of six harmonic modes of the functional data, 
F^(v,t) , as produced by program - MOD-5, one is faced with 
the retention of 18000 numbers which are required either 
in core -storage or in an on-line direct access device. In 
the IBM- 360/67, all real numbers require a 4 byte location, 
and a minimum problem would require 72k bytes (12 per node) 
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available during the running o£ the problem. Once the 
primary calculation o£ the density vector N(r,v,t) is done, 
data retention is required £or the total system N(r,v,t) 
which would be estimated as 4.0 (position points) (dimen- 
sions) (energy states) (time intervals) which can easily 
become a very large number when considering 40 or 50 points 
in a physical system, £or cases in two or three dimensions. 

The crux o£ the matter is that the one dimensional 
model will allo\v the utility o£ this approach to be demon- 
strated easily with relatively simple, easily compared 
Fourier series approximations in sine and cosine terms, 
rather than the £unctional £orms required £or solution in 
three dimensional cylindrical or spherical coordinate sys- 
tems £or sources not located at the principal axis symmetry. 

C. SOURCE FUNCTION APPROXIMATIONS 

Within the bounds o£ the one dimensional spatial model, 
the £ive source £unctions can be described by comparatively 
simple Fourier series approximations: 



£or a representation o£ the source geometry within the con- 
iines o£ the slab. 

The lollowing labeling conventions are used in delining 
the ternts ol the Fourier series representation: 

S^ source strength- total neutrons in the pulse. 




(A^ cos (k^x) + C„ sin (k x) ) 
n n n n 



n 



X 



o 



hall width ol the slab. 
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Series 

So urce Type Representation Fourier Coefficients 
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TABLE I. Fourier Expansion Coefficients. 
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D. 



ENERGY- VELOCITY -LETHARGY 



The sloiving down equation can be written in several 
forms to describe the variation of the neutron energy with 
time. For computational purposes, it is often convenient 
to describe the kinetic parameter in terms of the neutron 
velocity, while at other times, the kinetic energy of the 
neutron would provide a better description. There are three 
main terms that can be used interchangeably to describe 
the neutron during the slowing down process: velocity, 
energy and lethargy. Any expression in one of these param- 
eters can be transformed to the appropriate form in the 
others by the relationships: 

Energy-velocity E = m^v^/2 

Energy-lethargy - ln(E/E^) = u 

Velocity-lethargy - 21n(v/v^) = u. 

By these expressions, the parameters, E^ and v^ , de- 
scribe the neutron reference or source condition with the 
added convenience that as the neutron slows down, its leth- 
argy increases, while for the velocity and energy expres- 
sions, the slowing down process obviously results in the 
neutron going to a state with a smaller value. 

Simple relations between the derivatives of these ex- 
pressions will allow the slowing down equation to be easily 
revised with the desired terms: 

dE = m V dv 
n 
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du 



-2 dv 

V 



-dE 

E 



In terms of a neutron density function, N^Cv) and the 
differential relations, all three forms of the neutron den- 
sity are: 



N(v) 


dv = 


N(E)dE 


N(v) = 


m vN(E) 
n 


N(E) 


dE = 


N (u) du 


-E N(E) 


= N(u) 


N(v) 


du = 


N (u) du 


N(u) = 


^ N(v). 



(2m^E)'%(E) 
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IV. APPLICATIONS 



A. NEUTRON DENSITY 

The neutron density can be determined in the numerical 
model by the infinite series expression: 



such that the series representation converges to the analy- 
tical value as N increases to infinity. In the numerical 
method for this calculation, the series representation will 
be determined as a finite series, truncated to the first 3 
to 6 terms. In this form of the approximation, the error 
introduced is principally that of under estimating the true 
value if an even number of terms are used, while slightly 
over estimating with an odd number of terms. 

The truncation error of the finite series approximation 
can be estimated from the initial conditions of the system 
at time equal zero as the neutron density at that time must 
equal the initial source strength: 



If the total source strength (S^) is obtained by integrat- 
ing over- all velocities and the source volume, the source 
strength can be given by the series approximation: 



N 




1 



N(r ,v,0) = S(r,v,0) . 




V r 



V r n 
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With the use o£ a truncated series approximation, the 
error introduced can be calculated as: 

Error //I Fn(v,0) dr dv. 

1 

4 

This normalization error is used to obtain a first order 
adjustment factor, , by the simple approximation: 

A _ Error 

^0 ■■ S 

o 

and a normalization adjustment is determined as: 

Source Normalization = (1 + A„) S 

^ 0 o 

such that the truncated series will give a numerical value 
equal to the initially defined parameter -S^ . This now gives 
the series approximation as: 

N 

S(r,v,0) = XI (1 + AJ R^(r) F^(v,0) 

1 

where the artifice of the source strength normalization fac- 
tor (Aq) would be a succeedingly smaller correction as more 
terms are included in the finite series representation. 

This adjustment of the numerical value of the source 
strength has no effect on the velocity- time functions, 
F^(v,t), as this is a constant multiplier included in the 
spatial dependent expressions. 

In the discrete state, discrete time calculations of 
MOD-5, the probability that a neutron will be within a 
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specified velocity interval (v^) at the finite time step 
(t^), transforms the neutron density expression to the 
fo rm ; 



n 

and the total neutron population at any fixed time can be 
found by integration over the spatial coordinate and summa- 
tion over the velocity intervals; 



N(t^.) 



r i n 



B. NEUTRON FLUX 

The neutron flux, a scaler quantity to describe the 
net flow of neutrons per unit time and unit area, is cal- 
culated as the product of the speed of the neutron and the 
number density of neutrons having that speed: 

4>(r,v,t) = V N(r,v,t) 

and is approximated in this work as the discrete value: 
4)(r,v.,t.) = V. N(r,v.,t.)- 

v_> 1 » j 1 » 1 » y 

The total integrated flux can be analytically determined by 
an integration over all speeds, and in this work, is there- 
fore approximated by a summation over all discrete states 
of the density parameter N(r,v^,tj) to yield the result: 

4>(r,t.) = Z v^N(r,v^,tp = E Z 

i i n ^ 
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C. MEAN ENERGY OF NEUTRONS 

The spatial variation of the neutron mean energy 
provides valuable information on their migration and dif- 
fusion following the initial burst. This additional in- 
sight into the phenomena of "diffusion cooling" considers 
the general trend of the higher velocity components of the 
neutron population to leak out the boundries of the finite 
system. As such, the remaining neutrons, while still under- 
going interactions which remove them or lower their veloc- 
ity, have their total number decreased by this additional 
effect . 

In the discrete velocity state, discrete time model, 
the mean energy of the neutrons at a point r in the slab is 
determined as: 



In comparing the results of this mean energy determina- 
tion, two basic methods can be used: first, to study the 

variation in time of the mean energy at a fixed position, 
and secondly, to consider the ratio of the mean energy at 
all positions in the system, for all time periods, compared 
to a single reference position. Within this work, the mean 
energy is compared at each finite time period to the current 
value of the mean energy at the center of the slab, while 
the time variation of the mean energy at the reference point 
is followed in detail. 




E(r ,t^) 



1 




1 
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In the discrete time-velocity-space model, this is 



determined as: 



E(r,tj) 

E (r , t . ) 



^E(v.) N(r',v.,t.) 

J. 

ZN(r’,v.,t.) 

i 




i 



where r^ is the reference point, x=0.0, and r' is a position 
within the slab. 

D. DETECTOR RESPONSE 

Numerous methods [15] are currently utilized to deter- 
mine the neutron spectra and flux in fast reactor test as- 
semblies, for internal and external measurements. Fast 
neutron leakage is examined in time of flight tests to 
evaluate the time-varying spectra and neutron resonance 
absorption foils to study the internal spatial dependent 
spectra . 

For this analytical model, no single existing device 
appeared to be compatable to cover the continuous time and 
energy response as calculated by the discrete state method. 
An analytical detector is therefore constructed that would 
respond to the time- velocity variations of the neutron den- 
sity. For future applications of this model, a one dimen- 
sional version of an existing devices characteristics can 
be introduced. 
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The analytical detector is modeled as a proton recoil 
scintilation device capable of measuring the proton recoil 
of hydrogen-neutron collisions. The detector is considered 
to be filled with hydrogen at a molecular concentration as 
equivalent to an ideal gas at STP. , with physical dimen- 
sions of a width 2 about a midpoint 

The detector response (DR) is determined using the 
ABN-26 group cross section [16] to be: 



DR 



(r„,v.,t.) = f Z (v.) V. N(r',v.,t.) dr' 

D’ J 1-^ 1 — 



where the integral over the spatial coordinate dr' refers 
to the total detector volume. 

The total response of the detector both for the total 
integrated flux over all velocities and over all times can 
be determined in this discrete velocity- time model by nu- 
merical integration and summation: 



DRV(rp,t^) = Ys DR(rp,v^,t^) 
i 

DRT(ro,t.) = Y DRV(rp,t.) 

j 

where the final summation over all times from t=0 to t=tj 
provides the neutron fluence at the detector, 

E. MID- SIX COMPUTER PROGRAiM 

A general purpose computer program, MIL-SIX, was written 
in FORTRAN IV to process the time -energy (lethargy) response 
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date produced by MOD-5 to determine within this version o£ 
the transport approximation, the neutron density function- 
N(r,u,t) for the one dimensional infinite slab geometry. 
Several minor modifications were made in MOD-5 to provide 
the ability to calculate the time-energy response of a sys- 
tem with multiple harmonic modes and to provide the neutron 
probability density for each lethargy state at the same 
time steps. The only form of output system tested involved 
the use of punch card output from MOD-5, however, provisions 
have been made to employ nine-track magnetic tape and/or on- 
line disk storage capabilities (IBM-2314 units). Program 
MIL-SIX provides the general user the ability to accept any 
of the three forms output from MOD-5, for up to six (6) 
harmonic modes of data. The program can be modified easily 
to adjust for a higher order approximation using more than 
six modes with a simple change in the size of the storage 
arrays in the program. 

Program MIL-SIX is composed of the following subpro- 
grams and routines : ' 

Subroutine Function 

Main Principal control portion of the program; 

handle all logical decisions on which 
routines the program will perform 
INCONl Establishes all principal default param- 

eters for the program; reads all problem 
definition statements (punch cards) which 
determine the parameters that the program 
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will need to calculate the desired func- 



READIO 
READll 
E ST I MATE 1 

MODEL 1 



PHIl 



tions; provides the user with instruc- 
tions and information concerning actions 
the program will perform for the stated 
problem run. 

Processes all punch card output from 
program MOD- 5 . 

Processes all tape (9-track) from pro- 
gram MOD- 5 . 

Provides guidance information of the 
execution time to be expected for the 
running of a particular problem. 

Calculates the non-harmonic terms of the 
fourier space dependent functional ex- 
pression for the source geometry speci- 
fied for the particular problem; estimates 
the error introduced by the use of a 
finite number of terms in the fourier 
infinite series approximation; adjusts 
the values of the fourier non-harmonic 
coefficients to obtain a closer approxi- 
mation to the infinite series representa- 
tion; and provide a simple space dependent 
plot of the fourier expansion of the space 
dependent function on the on-line printer. 
Principal function sub-program of the pro- 
gram: to calculate the fourier space de- 



pendent harmonic terms in the calculation 



FOUTl (x. ) 



FOUT2 



of the terms of R (x) (n=the harmonic 
mode, x=the space point) for the one 
dimensional model. 

Principal determination subprogram to 
calculate the series of values for the 
function - N(x,u,t) = \ ‘^’^i ’ ^ j ^ 
for summed response due to all harmonic 
mode to give the neutron density for 
the single point - x^. 

This subprogram will calculate the two 

principal space dependent relationships: 

<j)(x, ,t.) and E(x, ,t. the integrated 
K 3 1 

(total) neutron flux and the mean energy 
of the neutrons at the point (x^) for a 
profile of 100 points across the slab 
for each time steps (t^), and will pro- 
vide a graphical plot of the functional 
forms : 

<}>(Xj^,tj) vs Xj^ 

E(Xk,t^) / E(0,tj) vs Xk 

via a on-line printing of a rough plot 
of the results or via the Cal-Comp Plotter 
for each set of discrete time step date 
(tj); also provides tabulated data on 
neutron flux, mean energy and mean ener- 
by ratio. 
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F0UT3 



F0UT4 



Detector response: this routine con- 

siders the simple model o£ a hydrogen 
recoil detector located inside the slab. 
This is not a true device, but one which 
determines the function: 

u n 

using the integrated position dependent 
flux (calculated at 5 points within the 
assumed detector position) ; a simple 
average is performed to obtain the flux 
at the detector, where the energy depen- 
dent microscopic cross section (Z^(u^) 
for STP hydrogen is based on the ABN-26 
group cross section tables. Starting 
with 4>(x, ,t.) - flux data, a numerical 

K j 

integration is performed to determine 
the integrated detector response. Out- 
put of the two functional forms of the 
detector response are given both by the 
on-line printer, and the cal-comp plotter. 
Provides a simple spectral plot of the 
data input from program MOD- 5, to show 
the distribution F (u.,t.) for each har- 
monic mode versus the lethargy state, 
showing the movement of the neutron dis- 
tribution to states of lower energy (higher 
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F0UT5 



lethargy) with time for each mode of 
harmonic buckling. 

Provides a more detailed spectral re- 
sponse plot of the MOD- 5 input data, by 
consideration of the fourier harmonic 
coefficient R f>^v) for the individual 
position point - x-^ weighting the har- 
monic data - to give the response spec- 
tra: 

a single plot for each harmonic mode of 

input data, and a total spectral response 

function N(x, ,u.,t.) = Y R (x.)? ("u.,t.), 
^ k’ i’ 1 ^ n^ k-^ n^ i’ 

summed over all modes and plotted versus 
lethargy state - u^, on the Cal-Comp 
plotter . 
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•FIGURE U FLOW CHART COMPUTER PROGRAM MIL-SIX 
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V. DISCUSSION OF RESULTS AND CONCLUSIONS 



A. PHYSICAL SYSTEMS 

The harmonic expansion method was examined for applica- 
tion to multiplying and non-multiplying systems, a beryllium 
slab and a modified form of the ZPR-3 assembly 6F. The 
physical parameters of the two systems are given in Table 
II. The nuclear characteristics were assumed to be: 



System 


Constituents 


Mass Density 
(gm/ cc) 


Total Modes 


Beryllium Slab 


Beryllium 


1.84 


6 


ZPR-3 (6F) 


Aluminum 


0.848 


3 




Iron* 


0.965 






U-235 


2.622 






U-238 


3.016 





(*Iron substituted for stainless steel.) 

A pulse of one 2.46 MeV neutron was fixed as a delta func- 
tion in time to study the decay of each harmonic component 
of the discrete state velocity- time response ^ ^ . Two 

versions of the slowing down process were examined: 

(1) follow each harmonic mode to the same final time 
(tj) as the fundamental mode, 

(2) follow decay of each harmonic mode to the time step 
at which' the probability density component for that mode was 
below 0;001 of the neutron remaining in the system. 

Both the beryllium and the assembly 6F systems v;ore 
studied to determine the variations in the neutron spectra 
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resulting from the inclusions of the higher harmonic modes. 
The persistance of the higher modes was observed in the 
early time history, but provided a negligible contribution 
toward the end of the decay time. 

B. DATA PRODUCTION 

The neutron probability density data for the lethargy- 
time response was calculated by MOD-5 for all harmonic modes 
using two limits on the systems considered. MOD-5 normally 
selects those discrete time steps to provide the probability 
density data under its internal selection criteria as though 
all problems were a fundamental mode type calculation. This 
freedom in selecting the output time steps was over ridden 
by modification of the existing control sequency in the basic 
MOD- 5 program. 

The fundamental mode calculation was done for a system 
to obtain the fixed time intervals and output time steps 
for the higher harmonic mode calculations. With the "lock- 
ing in" of the data output time steps, the higher harmonic 
modes were run to the same final discrete time as the funda- 
mental mode. Unfortunately, as the higher modes decay more 
rapidly, this resulted in the abnormal termination of several 
computer runs when the probability density dropped to values 
smaller than the IBM-360/67 is capable of processing, i.e., 
below 10 . As a result, all final computer runs were set 

to terminate when the probability density component for a 
mode dropped below an arbitrary limit of 0.001. 
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The time decay of the harmonic components is mainly 
affected by the leakage of neutrons from the system which 
is dependent on the buckling factor, which appears in 

the time dependent diffusion equation as a term: 
vD(v)B^N(x,v,t) . 

A summary of the MOD- 5 calculations is given in Table 
III for the two basic systems where all calculations were 
carried to the point where the probability that a neutron 
had leaked from the system, been absorbed or had dropped 
below the lowest energy was .999. In all cases, the higher 
harmonics were still providing a small contribution to the 
total population at times approaching 300ns, but the major 
contribution ^vas due mainly to the fundamental. 

C. ANALYSIS OF EXPANSION METHODS 

The adequacy of the truncated Fourier series to repre- 
sent the various source geometry conditions was evaluated 
by the comparison of the source strength parameter which was 
determined via the form: 

^ XT ‘ ■ 

/ o N 

^ (A^ Cos Kn X + Cn Sin Kn x) dx 

-X 1 

o 

N 

A = 1.0 - y C-) Sin Kn X - Cos Kn x . 

o A-. kn o kn o 

1 

The uncorrected source strength error was evaluated for 

the first 100 terms of the series expansion; which stabilized 
to almost constant- values after the first 30 to 40 terms. 
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System Mode Time Final Neutron Probability Components (Percentage) 

(ns) Slowing to Bottom Capture 

of Spectrum** Leakage non-fission Fission 
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TABLE III. Summary of MOD- 5 Computer Runs. 



Uncorrected Source Error (A ) 100 Term 

^ o-^ 





Source Type 


ZPR-3 (6F) 


Beryllium Slab 


1- 


Square Pulse 


0.0652 


’ 0.0265 


3- 


Exponential 


0.0231 


0.190 


5- 


Ramp Function 


0.0090 


0.0090 



The Fourier series expansions are defined such that 
with an infinite number of terms, the error should go to 
zero. However, calculations show that for the source func- 
tion approximations, this error goes to some non-zero value 
It should not be any surprise that these source terms have 
a non-zero error, since the Fourier approximation is based 
on the interval, where the distance ”6" represents 

the extrapolation distance correction for the boundry con- 
dition of forcing the neutron density to zero at that point 

In the source function testing, the physical limits of 
the slab are within a width 2x , rather than boundries of 
2(x^+6). This additional value of twice the extrapolation 
distance is responsible for apparent constant error. When 
the integration is done within the bounds of the slab, a 
loss occurs at the end points in the intervals' 

- fx +6) < x' < - X and +x > x > fx +6). 

^o-^- - o o- -^o-^ 

As the ratio 5/x^ is greater than zero for all real 
(non- inf inite systems) , this apparent difficulty in the ap- 
proximation will be almost negligible for systems with the 
slab width much greater than the extrapolation distance 
(6/x^<<1.0). 

Each of the three main source types (the square pulse, 
exponential and ramp functions) were studied to the limit o 
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a 100 mode harmonic expansion for the spatial dependent 

functions. The value of the delta correction fA ) are 

shown in Figures 5 and 6 for both systems. Ifliile all test 

calculations were done for a six mode expansion, the effect 

of the truncation to less than an infinite series approaches 

an asymtotic limit in each source type, the limits on the 

value of A are a measure of the 6/x ratio. In both of the 
o o 

systems studied, it can easily be seen that inclusion of 
more than 10 harmonic modes would not significantly change 
the spatial dependent values. 

D. COMPARISON OF SOURCE TYPES 

The development and testing of the time-space-energy' 
dependence for five source geometries was a task more for- 
midable than originally anticipated. Numerous programming 
difficulties arose toward the end of the project for the 
two cases of the off-axis sources, which resulted in their 
not being included in the final calculations. Both of these 
sources require modification of the spatial harmonic func- 
tions than those proposed previously. 

The three source conditions that are included in the 
detailed calculations were the symmetrical square pulse and 
the first collision source centered about the mid-point of 
the slab and the exterior wide beam or ramp function source. 

The existing version of the computer program, MIL-6, 
demonstrated the effectiveness of the three main source 
conditions, but the off-axis square pulse and first colli- 
sion geometry must be revised and retested prior to further 
analysis . 
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TOTAL MODES - RIRMONIC EXPANSION 




60 



TOTAL MODES - lilRi'IONIC EXPAfIFION 



E. NEUTRON SPECTRA 

The primary goal of this project was to demonstrate the 
spatial dependent response of the neutron spectral function: 



The contribution of each harmonic mode to the neutron 
spectra is determined by a Fourier space harmonic term, 
R^(x) , and the MOD-5 time - lethargy probability density for 
the n-th time mode, F^(v,t) . For each harmonic mode, the 
probability density function is subject to the basic con- 
dition : 



the correct source contribution. 

F. NEUTRON DENSITY AND NEUTRON FLUX 

The results of neutron density and neutron flux cal- 
culations are presented at a limited number of- points to 
be representative of the over all processes, however, in- 
clusion of all points in the slab would result in an inor- 
dinate number of graphs and tables of data. 

The neutron density, N(x, ,v.,t.)> plotted only for 

K 1 j 

the point Xj,=0.0 and a complete spectral plot is given in 
Figures 13 and 16, to be representative of the information 
possible to determine. All other calculations made by the 
program MIL-6 required the evaluation of the density at all 
points in the slab, and similar response curves were prepared. 




n 



F^(v,t) = 1.0 - (Absorption + Leakage) 
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The total neutron flux, 



“ Z ''i NCx^,v.,t.), 

i 

is presented for one location half way between the center 
and the surface of the slab, the same positions selected 
for the detector response calculations. Each plot is for 
a single system, either the Beryllium or assembly 6F, fol- 
lows the complete time decay of the total flux from one 
nanosecond to the time at which the neutron population had 
decayed to 0.11 of its initial value. In both cases, the 
first build up of the flux indicates the initial stages of 
the neutron motion \irithin the slab, to the point at 20ns 
at which the leakage from the free surface becomes the 
dominant effect. 

When the leakage terms are considered, the appearance 
of the decay curve now looks almost identical with the 
pseudo decay modes previously described. When looking at 
the neutron density curves, it can easily be noticed that 
it required several nanoseconds for the initial 2.46 MeV 
burst to develop into a smooth distribution. 

G. MEAN ENERGY AND MEAN ENERGY RATIO 

A convenient method of comparing the time and spatial 
effects on the neutron density is by comparing the average 
or mean energy of the neutrons at all points in the slab 
at the same time. A complete history of the neutron popu- 
lation can be concisely presented by observing the absolute 
variation of the mean energy at a single reference point in 
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the slab and then comparing the relative changes in the 
entire slab in terms of the reference point for each dis- 
crete time step. 

The basic reference point for all comparisons is the 
center of the slab, x=0.0, for all geometries and systems. 
The decay of the mean energy with time is shown for both 
beryllium (Figure 19) and assembly 6F (Figure 20) . Start- 
ing with the two reference cases, the spatial variation from 
the center to the surface of the slab starting at two nano- 
seconds after the pulse to the final decay of the funda- 
mental mode, (beryllium - Figures 21 to 30--assembly 6F 
Figures 31 to 36) . 

After the initial pulse, a wave-like shape is observed 
to develop across the slab in the direction of the free 
surface, as one would expect the higher energy neutrons to 
migrate out of the system fastest. This is shown most 
vividly in Figures 24 to 27. 

Near the end of the decay period, the last remaining 
neutrons in the system have almost a constant ratio of one, 
as the system tends toward an equilibrium condition. The 
leakage at the surface allows the final stages to show this 
condition (Figures 30 and 36) . 

H. DETECTOR RESPONSE 

This computational routine did not work as well as 
intended. Minor computational and programming errors were 
encountered very late in the project and have not been com- 
pletely resolved. These will require revision and further 
testing as follow on work. 
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Comparisons of the flux calculations and the mean en- 
ergy determinations referred to in Sections F and G were 
compared to the "proton recoil" detector (Figures 37 and 38) 
which follow the same basic trends for the time dependent 
response, but the numerical differences in the time inte- 
grated response will require improvement of the space in- 
tegral and time summation numerical methods. 
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VI . RECOMMENDATIONS 



A. FOLLOW ON WORK 

The program MIL-6 was prepared to be a general compu- 
tational tool using the harmonic mode expansion method for 
the one dimensional systems. Principal follow on work is 
recommended to correct the minor computational difficul- 
ties encountered, v\rrite and test data transfer routines 
for both MOD-5 and MIL-6 and review and edit the program- 
ming for a more optimized execution. 

Principal areas to have future work are: 

(1) re-evaluate the Fourier space functions for the 

two off-axis source geometries, , 

(2) revise and test the numerical integration tech- 
niques of the detector response program for the summation 
over time, and consider methods to improve the space inte- 
gration over the detector volume, 

(3) consider a detector response function that would 
be more realistic than the STP hydrogen case now included 
in the program, 

(4) complete and test the input/output options using 
nine-track magnetic tape and the IBM-2314 magnetic disc 
systems for data transfer between the programs MOD-5 and 
MIL-6, 

(5) revise the program MIL-6 to handle more than six 
harmonic modes by use of a smaller state structure of 70 
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to 100 states rather than the 150 states the program is 
currently written to process, 

(6) provide more realistic source geometries for the 
one dimensional case than the five simple ones currently 
in the program, 

(7) review the program for reduction in the core 
storage space requirements and to improve the speed of 
execution through program optimization. 

B. IMPROVEMENTS TO PROGRAM MIL- 6 

The original program MIL-6 was written to deal with 
monoenergetic neutron sources and delta functions in time. 
This may require significant modification of the computa- 
tional methods to considering the fission spectrum sources 
and sources that have a finite width time duration. These 
situations had not been tested in the MIL-6 work and would 
provide a more realistic model for fast neutron experiments 

A second item that appears worthy of consideration in 
the future will assist in modeling the wave-like motion of 
the neutron pulse from the finite volume square pulse and 
the exterior wide beam ramp function. This could be con- 
sidered as a limiting correction factor to couple the time- 
space and energy dependence during the initial response to 
a neutron pulse. 

C. THREE DIMENSIONAL SYSTEMS 

The harmonic expansion method has been demonstrated to 
work for the one dimensional case, and definite action to 



66 



begin study of three dimensional systems is warranted. A 
recommended sequence of priorities would be to expand the 
MIL-6 work to three dimensions in the cartesian coordi- 
nates first, cylindrical geometry second and spherical 
geometry coordinate systems last. 
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APPENDIX A: NUMERICAL RESULTS 



A. SPECTRAL RESPONSE 
1 . Beryllium Slab 

a. Neutron density harmonic component N^(x,v^,tj) 
for harmonic modes 1 to 6 at point x=0.0. (Figures 7 to 
12 .) 



b. Neutron density: N(x,v^,tj): lethargy - 

time neutron population density versus lethargy for a 6 
mode expansion at point x=0.0. (Figure 13.) 

2 . ZPR-5 Assembly 6F 

a. Neutron density harmonic component N^(x,Vj|^ ,tj ) 
for harmonic modes - 1 and 2 at point x=0.0. (Figures 14 
and 15.) 



b. Neutron density N(x,v^,t^) - lethargy - time 
neutron population density versus lethargy - for a sum of 
3 harmonic modes at x=0.0. (Figure 16.) 



B. MEAN ENERGY AND MEAN ENERGY RATIO 
1 . Beryllium 

a. Mean energy verus time at x=0.0. (Figure 19.) 

b. Plot of mean energy ratio versus position for 
8 time periods. (Figures 21 to 30) 



2. ZPR-3 (6F) 

a. Mean energy versus time at x=0.0. (Figure 20.) 

b. Plot of mean energy ratio versus position for 
6 time periods. (Figures 31 to 36.) 
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C. DETECTOR RESPONSE 

1. Beryllium - Figure 37 

2. ZPR-3 - Figure 38 
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FIGURE 26 
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FIGURE 35 
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NTMIN=NTVAL-1 

0X1=( DUMMY( NTVAL )+DUMMY( NTMIN) )^ 0.50 
DX2 = TIM( NIVAL )-TIM(I\ITMlN) 

DX3 = DX1^<0X2 
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